library(fixest)
library(modelsummary)

options(scipen = 999) ##scientific notation (no)

data_2021 <- readRDS("data_2021.rds")  
data_2024 <- readRDS("data_2024.rds")  

########################################################################################################################################################################################
########################################################################################################################################################################################
########################################################################################################################################################################################

model1_2024 <- feols(Q18_2_clean ~ 1|p_ottoman ~ route_density_2 + wheat_rainfed + istanbul_distance + vienna_distance, data = data_2024, vcov = "conley")
model2_2024 <- feols(Q18_2_clean ~ D1 + D2 + D3 + D4 + D5A + D6|p_ottoman ~ route_density_2 + wheat_rainfed + istanbul_distance + vienna_distance, data = data_2024, vcov = "conley")
model3_2024 <- feols(Q18_2_clean ~ D1 + D2 + D3 + D4 + D5A + D6 + Q2 + Q5|p_ottoman ~ route_density_2 + wheat_rainfed + istanbul_distance + vienna_distance, data = data_2024, vcov = "conley")
model4_2024 <- feols(Q19_2_clean ~ 1|p_ottoman ~ route_density_2 + wheat_rainfed + istanbul_distance + vienna_distance, data = data_2024, vcov = "conley")
model5_2024 <- feols(Q19_2_clean ~ D1 + D2 + D3 + D4 + D5A + D6|p_ottoman ~ route_density_2 + wheat_rainfed + istanbul_distance + vienna_distance, data = data_2024, vcov = "conley")
model6_2024 <- feols(Q19_2_clean ~ D1 + D2 + D3 + D4 + D5A + D6 + Q2 + Q5|p_ottoman ~ route_density_2 + wheat_rainfed + istanbul_distance + vienna_distance, data = data_2024, vcov = "conley")

# Summary for 2024 data
msummary(list(model1_2024, model2_2024, model3_2024, model4_2024, model5_2024, model6_2024),
         stars =c('.' = .1,'*' = .5, '**' = .01, '***' = .001))

########################################################################################################################################################################################
########################################################################################################################################################################################
########################################################################################################################################################################################

model1_2021 <- feols(Q18_2 ~1|p_ottoman ~ route_density_2 + wheat_rainfed + istanbul_distance + vienna_distance, data = data_2021, vcov = "conley")
model2_2021 <- feols(Q18_2 ~  + D1 + D2 + D3 + D4 + D5A + D6|p_ottoman ~ route_density_2 + wheat_rainfed + istanbul_distance + vienna_distance, data = data_2021, vcov = "conley")
model3_2021 <- feols(Q18_2 ~  + D1 + D2 + D3 + D4 + D5A + D6 + Q2 + Q5|p_ottoman ~ route_density_2 + wheat_rainfed + istanbul_distance + vienna_distance, data = data_2021, vcov = "conley")
model4_2021 <- feols(Q19_2 ~ 1|p_ottoman ~ route_density_2 + wheat_rainfed + istanbul_distance + vienna_distance, data = data_2021, vcov = "conley")
model5_2021 <- feols(Q19_2 ~  + D1 + D2 + D3 + D4 + D5A + D6|p_ottoman ~ route_density_2 + wheat_rainfed + istanbul_distance + vienna_distance, data = data_2021, vcov = "conley")
model6_2021 <- feols(Q19_2 ~  + D1 + D2 + D3 + D4 + D5A + D6 + Q2 + Q5|p_ottoman ~ route_density_2 + wheat_rainfed + istanbul_distance + vienna_distance, data = data_2021, vcov = "conley")

# Summary for 2021 data
msummary(list(model1_2021, model2_2021, model3_2021, model4_2021, model5_2021, model6_2021),
         stars =c('.' = .1,'*' = .5, '**' = .01, '***' = .001))

